Modelling of axonal locasation
Raphaelle Luisier, Idiap Research Institute
February 2023
1 Load the data
3' UTR isoform quantification from 8 samples originating from rat neuronal sympathetic cell cultures where the cell body is treated with low does of NGF and distal axons are either exposed to NGF or NT3. There are a total of 8 samples: cell body versus distal axons in NGF versus NT3 culture condition across 2 technical replicates.
We first need to select the reliably expressed genes and transcripts i.e. those which are reliably expressed in the cell body
par(mfrow=c(2,4))
coloi <- c("NGF.axon.1.raw.corrected","NGF.axon.2.raw.corrected","NGF.cb.1.raw.corrected","NGF.cb.2.raw.corrected","NT3.axon.1.raw.corrected","NT3.axon.2.raw.corrected","NT3.cb.1.raw.corrected","NT3.cb.2.raw.corrected")
mydat <- anno_tot[which(anno_tot$is.conservative),match(coloi,colnames(anno_tot))]
mydatall <- anno_tot[,match(coloi,colnames(anno_tot))]
for(i in c(1:ncol(mydat))){
mydat[,i] <- log2(mydat[,i]+1)
}
lims <- apply(mydat,2,function(Z)2^SelectExpressed(dat=Z,frac.bg=0.95,frac.fg=0.1))
tempsel <- apply(mydatall,2,function(Z)return(Z>=mean(lims[1,c(1,2)])))
soft.sel <- cbind(tempsel[,1]&tempsel[,2],tempsel[,3]&tempsel[,4],tempsel[,5]&tempsel[,6],tempsel[,7]&tempsel[,8])
tempsel <- do.call(what=cbind,lapply(c(1:8),function(Z)return(mydatall[,Z]>=lims[2,Z])))
hard.sel <- cbind(tempsel[,1]|tempsel[,2],tempsel[,3]|tempsel[,4],tempsel[,5]|tempsel[,6],tempsel[,7]|tempsel[,8])
final.sel <- cbind(soft.sel[,1]|hard.sel[,1],soft.sel[,2]|hard.sel[,2],soft.sel[,3]|hard.sel[,3],soft.sel[,4]|hard.sel[,4])
anno_tot$NGF.axon.is.expressed.iso <- final.sel[,1]
anno_tot$NGF.cb.is.expressed.iso <- final.sel[,2]
anno_tot$NT3.axon.is.expressed.iso <- final.sel[,3]
anno_tot$NT3.cb.is.expressed.iso <- final.sel[,4]
#Select only those reliably expressed in CB of both NGF and NT3
coloi <- grep(colnames(anno_tot),pattern=".raw.corrected")
seloi <- anno_tot$NT3.cb.is.expressed.iso&anno_tot$NGF.cb.is.expressed.iso#If we want to study differential transport, we need to make sure that things are expressed at least in both NGF and NT3 CB
#seloi <- c(1:nrow(anno_tot))
data.con.cb <- anno_tot[seloi,coloi]
row_info <- anno_tot[seloi,]
anno_tot <- anno_tot[seloi,]
mydatall <- log2(1+data.con.cb)
myUTR <- myUTR[seloi,]#B. Compute average expression in each compartment
avg.ngf.axons <- apply(mydatall[,c("NGF.axon.1.raw.corrected","NGF.axon.2.raw.corrected")],1,mean)
avg.nt3.axons <- apply(mydatall[,c("NT3.axon.1.raw.corrected","NT3.axon.2.raw.corrected")],1,mean)
avg.ngf.cb <- apply(mydatall[,c("NGF.cb.1.raw.corrected","NGF.cb.2.raw.corrected")],1,mean)
avg.nt3.cb <- apply(mydatall[,c("NT3.cb.1.raw.corrected","NT3.cb.2.raw.corrected")],1,mean)
avg.cb <- apply(mydatall[,c("NT3.cb.1.raw.corrected","NT3.cb.2.raw.corrected","NGF.cb.1.raw.corrected","NGF.cb.2.raw.corrected")],1,mean)
avg.ax <- apply(mydatall[,c("NGF.axon.1.raw.corrected","NGF.axon.2.raw.corrected","NT3.axon.1.raw.corrected","NT3.axon.2.raw.corrected")],1,mean)
avg.ngf.axons.raw <- apply(data.con.cb[,c("NGF.axon.1.raw.corrected","NGF.axon.2.raw.corrected")],1,mean)
avg.nt3.axons.raw <- apply(data.con.cb[,c("NT3.axon.1.raw.corrected","NT3.axon.2.raw.corrected")],1,mean)
avg.ngf.cb.raw <- apply(data.con.cb[,c("NGF.cb.1.raw.corrected","NGF.cb.2.raw.corrected")],1,mean)
avg.nt3.cb.raw <- apply(data.con.cb[,c("NT3.cb.1.raw.corrected","NT3.cb.2.raw.corrected")],1,mean)2 Differential transport using classic abundance ratios
#A. Compute relative log ratio
RUD.ngf <- avg.ngf.axons-avg.ngf.cb
RUD.nt3 <- avg.nt3.axons-avg.nt3.cb
diff.RUD <- RUD.ngf-RUD.nt3
#B. Compute relative usage (PUD has larger dynamic range than SUD; with SUD lowly expressed genes in CB have higher chances of getting high score so better PUD which detects changes across the entire spectrum of CB expression level)
PUD.ngf <- avg.ngf.axons/(avg.ngf.cb+avg.ngf.axons)
PUD.nt3 <- avg.nt3.axons/(avg.nt3.cb+avg.nt3.axons)
PUD.nt3[(avg.nt3.cb+avg.nt3.axons)==0] <-0
PUD.ngf[(avg.ngf.cb+avg.ngf.axons)==0] <-0
diff.PUD <- PUD.ngf-PUD.nt3
SUD.ngf <- avg.ngf.axons.raw/(avg.ngf.cb.raw+avg.ngf.axons.raw)
SUD.nt3 <- avg.nt3.axons.raw/(avg.nt3.cb.raw+avg.nt3.axons.raw)
SUD.nt3[(avg.nt3.cb.raw+avg.nt3.axons.raw)==0] <-0
SUD.ngf[(avg.ngf.cb.raw+avg.ngf.axons.raw)==0] <-0
diff.SUD <- SUD.ngf-SUD.nt3
LUD.ngf <- avg.ngf.axons.raw/avg.ngf.cb.raw
LUD.nt3 <- avg.nt3.axons.raw/avg.nt3.cb.raw
LUD.nt3[avg.nt3.cb.raw==0] <-0
LUD.ngf[avg.ngf.cb.raw==0] <-0
diff.LUD <- LUD.ngf/SUD.nt3
#C. Compute Fisher count P-value --> P-value on differential transport between NGF and NT3 while previous analysis is about transport efficiency
mytempdat <- data.con.cb
sum.dat <- t(apply(mytempdat,1,function(Z)return(tapply(Z,INDEX=c("NGF.axon","NGF.axon","NGF.cb","NGF.cb","NT3.axon","NT3.axon","NT3.cb","NT3.cb"),FUN=sum))))
p.val.diff <- apply(sum.dat,1,function(Z)return(fisher.test(round(cbind(Z[c(1,2)],Z[c(3,4)])))$p.value))
FDR.diff <- p.adjust(p.val.diff,method="fdr")
#D. Prepare output
PUD <- data.frame(PUD.ngf=PUD.ngf,PUD.nt3=PUD.nt3)
RUD <- data.frame(RUD.ngf=RUD.ngf,RUD.nt3=RUD.nt3)
SUD <- data.frame(SUD.ngf=SUD.ngf,SUD.nt3=SUD.nt3)
dPUD.rel <- PUD$PUD.ngf*(avg.ngf.cb/max(avg.ngf.cb))-PUD$PUD.nt3*(avg.nt3.cb/max(avg.nt3.cb))
dSUD.rel <- SUD$SUD.ngf*(avg.ngf.cb/max(avg.ngf.cb))-SUD$SUD.nt3*(avg.nt3.cb/max(avg.nt3.cb))
myOut <- data.frame(anno_tot[,c("txID","geneSymbol","uniqueID")],FDR=FDR.diff,dRUD=diff.RUD,dPUD=diff.PUD,dSUD=diff.SUD,dPUD.rel=dPUD.rel,dSUD.rel=dSUD.rel,PUD,SUD,
PUD.ngf.rel=PUD$PUD.ngf*(avg.ngf.cb/max(avg.ngf.cb)),
SUD.ngf.rel=SUD$SUD.ngf*(avg.ngf.cb/max(avg.ngf.cb)),
PUD.nt3.rel=PUD$PUD.nt3*(avg.nt3.cb/max(avg.nt3.cb)),
SUD.nt3.rel=SUD$SUD.nt3*(avg.nt3.cb/max(avg.nt3.cb)),
dRUD=diff.RUD,RUD,data.con.cb)
#G. Extract differentially transported txID based on RUD and PUD
#Test for overlap with Repetitive regions -- test whether this is required
myRpM <- import.gff("./zenodo/rmsk_rn5_ucsc.sorted.gff",format="gff")
myRpM <- reduce(myRpM,min.gapwidth=20)
POS <- as.character(strand(myUTR))=="+"
myUTR.focus <- myUTR
start(myUTR.focus)[POS] <- end(myUTR)[POS]-499
end(myUTR.focus)[!POS] <- start(myUTR)[!POS]+499
gOver <- findOverlaps(query=myUTR.focus,subject=myRpM,ignore.strand=FALSE)
idxU <- queryHits(gOver)
idxG <- subjectHits(gOver)
new.start <- apply(cbind(start(myUTR.focus)[idxU],start(myRpM)[idxG]),1,max)
new.end <- apply(cbind(end(myUTR.focus)[idxU],end(myRpM)[idxG]),1,min)
my.overlapping.utr <- myUTR.focus[idxU,]
start(my.overlapping.utr) <- new.start
end(my.overlapping.utr) <- new.end
mywidth <- tapply(width(my.overlapping.utr),INDEX=factor(as.character(my.overlapping.utr$ID)),FUN=sum)
oi.rpm <- names(mywidth)[mywidth>=300]
#save(list="myOut",file="./zenodo/transport/diff_transport_ratios/Differential_abundance_ratios.RData")
sel.A.up <- myOut$dRUD>1
sel.A.do <- myOut$dRUD<(-1)
sel.B.up <- myOut$dPUD.rel>0.10
sel.B.do <- myOut$dPUD.rel<(-0.10)
sel.C <- myOut$FDR<0.01
sel.D <- !(myOut$uniqueID%in%oi.rpm)
sel.more.NGF <- sel.A.up&sel.B.up&sel.D&sel.C&anno_tot$NGF.axon.is.expressed.iso&data.con.cb[,1]>20&data.con.cb[,2]>20
sel.more.NT3 <- sel.A.do&sel.B.do&sel.D&sel.C&anno_tot$NT3.axon.is.expressed.iso&(data.con.cb[,5]>20|data.con.cb[,6]>20)#More laxist on NT3 since it is not as covered as the NGF
txID.more.NGF <- unique(as.character(anno_tot$txID)[sel.more.NGF])#643
txID.more.NT3 <- unique(as.character(anno_tot$txID)[sel.more.NT3])#322
#write.csv(myOut[sel.more.NGF,],"./zenodo/transport/diff_transport_ratios/more_in_NGF.csv")
#write.csv(myOut[sel.more.NT3,],"./zenodo/transport/diff_transport_ratios/more_in_NT3.csv")more_NGF <- read.csv("./zenodo/transport/diff_transport_ratios/more_in_NGF.csv")
more_NT3 <- read.csv("./zenodo/transport/diff_transport_ratios/more_in_NT3.csv")
load("./zenodo/transport/diff_transport_ratios/Differential_abundance_ratios.RData")3 Development of the axonal localisation score
3.1 Selection of linear model that best fit the data
We next want to test whether the read count in the axonal compartment depends on the read count in the cell body as well as on transcript length (log10). We do this by fitting a linear model of iteratively increasing order for NGF and NT3 and then extract AIC.
Below is the result for the NGF condition:
# C. Modelling axonal transport in function of 1) read counts in CB and 2) tx length
# With the linear model --using both size and expression in CB -- the aim is to identify those which are not behaving in the same way as most other species of the same size
# Linear model to measure the extend of deviation from the linear model
dat.ngf <- data.frame(txL=log10(anno_tot$txLength[match(names(avg.ngf.axons),anno_tot$uniqueID)]),
cb=avg.ngf.cb,
axons=avg.ngf.axons)
dat.nt3 <- data.frame(txL=log10(anno_tot$txLength[match(names(avg.nt3.axons),anno_tot$uniqueID)]),
cb=avg.nt3.cb,
axons=avg.nt3.axons)
dat.ngf <- data.frame(txL=log10(anno_tot$txLength),
cb=avg.ngf.cb,
axons=avg.ngf.axons)
dat.nt3 <- data.frame(txL=log10(anno_tot$txLength),
cb=avg.nt3.cb,
axons=avg.nt3.axons)
polyfit <- function(i) x <- AIC(lm(axons~poly(cb,i)+txL,dat=dat.ngf))
fit1.1 <- lm(axons~cb,dat=dat.ngf)
fit1.2 <- lm(axons~cb+txL,dat=dat.ngf)
fit2.1 <- lm(axons~poly(cb,2),dat=dat.ngf)
fit2.2 <- lm(axons~poly(cb,2)+txL,dat=dat.ngf)
fit3.1 <- lm(axons~poly(cb,3),dat=dat.ngf)
fit3.2 <- lm(axons~poly(cb,3)+txL,dat=dat.ngf)
fit4.1 <- lm(axons~poly(cb,4),dat=dat.ngf)
fit4.2 <- lm(axons~poly(cb,4)+txL,dat=dat.ngf)#selected one
fit5.1 <- lm(axons~poly(cb,5),dat=dat.ngf)
fit5.2 <- lm(axons~poly(cb,5)+txL,dat=dat.ngf)
print(AIC(fit1.1,fit1.2,fit2.1,fit2.2,fit3.1,fit3.2,fit4.1,fit4.2,fit5.1,fit5.2))## df AIC
## fit1.1 3 111813.1
## fit1.2 4 110811.2
## fit2.1 4 109015.4
## fit2.2 5 108049.9
## fit3.1 5 108812.1
## fit3.2 6 107837.8
## fit4.1 6 108796.7
## fit4.2 7 107825.3
## fit5.1 7 108798.5
## fit5.2 8 107827.0print(anova(fit4.2,fit5.2,test="Chisq"))## Analysis of Variance Table
##
## Model 1: axons ~ poly(cb, 4) + txL
## Model 2: axons ~ poly(cb, 5) + txL
## Res.Df RSS Df Sum of Sq Pr(>Chi)
## 1 27654 79825
## 2 27653 79825 1 0.85188 0.587print(anova(fit4.2,fit4.1,test="Chisq"))## Analysis of Variance Table
##
## Model 1: axons ~ poly(cb, 4) + txL
## Model 2: axons ~ poly(cb, 4)
## Res.Df RSS Df Sum of Sq Pr(>Chi)
## 1 27654 79825
## 2 27655 82685 -1 -2859.5 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Below is the result for the NT3 condition.
fit1.1 <- lm(axons~cb,dat=dat.nt3)
fit1.2 <- lm(axons~cb+txL,dat=dat.nt3)
fit2.1 <- lm(axons~poly(cb,2),dat=dat.nt3)
fit2.2 <- lm(axons~poly(cb,2)+txL,dat=dat.nt3)
fit3.1 <- lm(axons~poly(cb,3),dat=dat.nt3)
fit3.2 <- lm(axons~poly(cb,3)+txL,dat=dat.nt3)
fit4.1 <- lm(axons~poly(cb,4),dat=dat.nt3)
fit4.2 <- lm(axons~poly(cb,4)+txL,dat=dat.nt3)#selected one
fit5.1 <- lm(axons~poly(cb,5),dat=dat.nt3)
fit5.2 <- lm(axons~poly(cb,5)+txL,dat=dat.nt3)
print(AIC(fit1.1,fit1.2,fit2.1,fit2.2,fit3.1,fit3.2,fit4.1,fit4.2,fit5.1,fit5.2))## df AIC
## fit1.1 3 104426.29
## fit1.2 4 103368.63
## fit2.1 4 100569.12
## fit2.2 5 99526.10
## fit3.1 5 100443.77
## fit3.2 6 99387.70
## fit4.1 6 100412.76
## fit4.2 7 99362.37
## fit5.1 7 100414.63
## fit5.2 8 99364.19print(anova(fit4.2,fit5.2,test="Chisq"))## Analysis of Variance Table
##
## Model 1: axons ~ poly(cb, 4) + txL
## Model 2: axons ~ poly(cb, 5) + txL
## Res.Df RSS Df Sum of Sq Pr(>Chi)
## 1 27654 58785
## 2 27653 58784 1 0.37677 0.6738print(anova(fit4.2,fit4.1,test="Chisq"))## Analysis of Variance Table
##
## Model 1: axons ~ poly(cb, 4) + txL
## Model 2: axons ~ poly(cb, 4)
## Res.Df RSS Df Sum of Sq Pr(>Chi)
## 1 27654 58785
## 2 27655 61064 -1 -2279.7 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1myfits <- list(fit1.1,fit1.2,fit2.1,fit2.2,fit3.1,fit3.2,fit4.1,fit4.2,fit5.1,fit5.2)
#par(mfrow=c(5,2))
#lapply(myfits,function(Z){
# pdt <- predict(Z)
# plot(avg.ngf.cb,avg.ngf.axons,pch=19,col=rgb(0,0,0,0.1),cex=0.2,main=,cex.main=0.7)
# points(avg.ngf.cb,pdt,pch=19,col=rgb(1,0,0,0.1),cex=0.2)
#})The selected model is polynom of degree 4 with transcript length. Here we check that the predictions (red dots) are close to the observed read counts for the selected model.
par(mfrow=c(1,1))
Z=myfits[[8]]
pdt <- predict(Z)
plot(avg.ngf.cb,avg.ngf.axons,pch=19,col=rgb(0,0,0,0.1),cex=0.2,cex.main=0.7,las=1,frame=FALSE,main="NGF")
points(avg.ngf.cb,pdt,pch=19,col=rgb(1,0,0,0.1),cex=0.2)
Interestingly, transcripts that are highly expressed in cell bodies are not necessarily enriched in axons or dendrites, indicating that at least a subset of RNAs do not reach the peripheral compartments by 'passive' transport but are sorted and delivered with an selective mechanism.
3.2 Compute the localisation score
We next use use a Bayesian approach to get a localisation score with positive values indicating over-transport and negative values indicating under-transport. We first group 3' UTR isoforms of similar expression level in the CB and of similar transcript length together. Assuming that isoforms within a group should exhibit similar expression in the cell body, of mean and standard-deviation derived from the linear model learned in the previous step. For each 3' UTR isoform, we calculate the likelihood of that the observed axonal read count are higher or lower than the predicted read counts given the prior distribution.
#Extract the probability of each data-point being generated to identify those which are under- or over-transported
#A. Fit polynomial of degree 4 to each data-set
fit4.2.ngf <- lm(axons~poly(cb,4)+txL,dat=dat.ngf)#selected one
fit4.2.nt3 <- lm(axons~poly(cb,4)+txL,dat=dat.nt3)#selected one
#B. Model prediction
pdt.ngf <- predict(fit4.2.ngf)
pdt.nt3 <- predict(fit4.2.nt3)
plot(pdt.ngf,pdt.nt3,pch=19,col=rgb(0,0,0,0.2),cex=0.2)
RUD.ngf <- avg.ngf.axons-avg.ngf.cb
RUD.nt3 <- avg.nt3.axons-avg.nt3.cb
out.ngf <- data.frame(dat.ngf,prediction=pdt.ngf,RUD=RUD.ngf,anno_tot[match(names(avg.ngf.axons),anno_tot$uniqueID),],avg.ngf.cb,avg.ngf.axons)
out.nt3 <- data.frame(dat.nt3,prediction=pdt.nt3,RUD=RUD.nt3,anno_tot[match(names(avg.nt3.axons),anno_tot$uniqueID),],avg.nt3.cb,avg.nt3.axons)
#D. Create intervals of expression according to expression in CB and transcript length
CreateBins <- function(myout){
bins.cb <- cut(myout$cb,breaks=seq(from=3,to=20,by=1),iclude.lowest=T)
bins.txL <- cut(myout$txL,breaks=seq(from=2,to=4.5,by=0.25),iclude.lowest=T)
bins.both <- paste(bins.cb,bins.txL,sep=".")
temp <- data.frame(table(bins.both))
tokeep <- temp$bins.both[temp$Freq>=5]
#tokeep <- temp$bins.both[temp$Freq>=1]
bins.both[!bins.both%in%tokeep]<-NA
return(factor(as.character(bins.both)))
}
bins.both.ngf <- CreateBins(myout=out.ngf)
bins.both.nt3 <- CreateBins(myout=out.nt3)
#E.Get the probability of the observed read count in the axons larger or smaller than the observed read counts in axons given their observed read counts in the CB and their transcript length
GetProbs <- function(mybins=bins.both.ngf,myfit=fit4.2.ngf,out=out.ngf){
myOut <- do.call(lapply(c(1:length(levels(mybins))),function(IX){
print(IX)
myprediction <- compute.prediction(IX,myfit=myfit,bins=mybins)
subdat <- out[which(mybins==levels(mybins)[IX]),]
if(nrow(subdat)>=10){
#Fit prior normal distribution of the read counts in the axons with prior mean and sd from prediction given the intervals
fit.norm <- fitdist(subdat$axons, "norm",start=list(mean=mean(myprediction),sd=sd(myprediction)))
#Sample from the prior
my.rdm.norm <- rnorm(n=10^4, mean=fit.norm$estimate["mean"], sd=fit.norm$estimate["sd"])
#Compute the posterior prob. i.e. empirical P-value
temp <- do.call(what=rbind,lapply(subdat$axons,function(Z)return(c(sum(my.rdm.norm<=Z)/10^4,sum(my.rdm.norm>=Z)/10^4))))
colnames(temp)<- c("prob(N).smaller","prob(N).bigger")
temp <-data.frame(temp,subdat)
return(temp)
}
else{
temp <- matrix(NA,ncol=2,nrow=nrow(subdat))
colnames(temp)<- c("prob(N).smaller","prob(N).bigger")
temp <-data.frame(temp,subdat)
return(temp)
}
}),what=rbind)
p.value.norm <- apply(cbind(-log10(myOut$prob.N..smaller+min(myOut$prob.N..smaller[myOut$prob.N..smaller>0],na.rm=TRUE)),
-log10(myOut$prob.N..bigger+min(myOut$prob.N..bigger[myOut$prob.N..bigger>0],na.rm=TRUE))),1,function(Z){
M<- max(Z)
if(is.na(M))return(NA)
if(which(Z==M)[1]==1)return(-M)
else
return(M)
})
myOut <- data.frame(p.value.norm,myOut)
return(myOut)
}
myOut.ngf <- GetProbs(mybins=bins.both.ngf,myfit=fit4.2.ngf,out=out.ngf)
myOut.nt3 <- GetProbs(mybins=bins.both.nt3,myfit=fit4.2.nt3,out=out.nt3)
ix <- match(rownames(myOut.nt3),rownames(myOut.ngf))
ix <- ix[!is.na(ix)]
myOut.ngf <- myOut.ngf[ix,]
ix <- match(rownames(myOut.ngf),rownames(myOut.nt3))
ix <- ix[!is.na(ix)]
myOut.nt3 <- myOut.nt3[ix,]
transport <- data.frame(uniqueIDs=rownames(myOut.ngf),
GS=as.character(anno_tot$geneSymbol[match(rownames(myOut.ngf),anno_tot$uniqueID)]),
txID=as.character(myOut.ngf$txID),
diff.transport=myOut.ngf$p.value.norm-myOut.nt3$p.value.norm,
diff.RUD=myOut.ngf$RUD-myOut.nt3$RUD,
transport.ngf=myOut.ngf$p.value.norm,
transport.nt3=myOut.nt3$p.value.norm,
RUD.ngf=myOut.ngf$RUD,
RUD.nt3=myOut.nt3$RUD,
cb.ngf=myOut.ngf$cb,
cb.nt3=myOut.nt3$cb,
axons.ngf=myOut.ngf$axons,
axons.nt3=myOut.nt3$axons,
anno_tot$NGF.cb.1.raw.corrected[match(rownames(myOut.ngf),anno_tot$uniqueID)],
anno_tot$NGF.cb.2.raw.corrected[match(rownames(myOut.ngf),anno_tot$uniqueID)],
anno_tot$NGF.axon.1.raw.corrected[match(rownames(myOut.ngf),anno_tot$uniqueID)],
anno_tot$NGF.axon.2.raw.corrected[match(rownames(myOut.ngf),anno_tot$uniqueID)],
anno_tot$NT3.cb.1.raw.corrected[match(rownames(myOut.ngf),anno_tot$uniqueID)],
anno_tot$NT3.cb.2.raw.corrected[match(rownames(myOut.ngf),anno_tot$uniqueID)],
anno_tot$NT3.axon.1.raw.corrected[match(rownames(myOut.ngf),anno_tot$uniqueID)],
anno_tot$NT3.axon.2.raw.corrected[match(rownames(myOut.ngf),anno_tot$uniqueID)])
transport$txLength <- anno_tot$txLength[match(rownames(myOut.ngf),anno_tot$uniqueID)]
transport <- transport[!is.na(transport$transport.ngf),]
mygr <- myUTR[match(transport$uniqueIDs,myUTR$ID),]
transport$coord <- paste(seqnames(mygr),":",start(mygr),"-",end(mygr),sep="")
tempax <- do.call(what=rbind,args=lapply(c(1:nrow(transport)),function(ix){
sel=transport$txID==transport$txID[ix]&transport$wL>transport$wL[ix]
return(c(ifelse(sum(sel)>0,max(transport$transport.ngf[sel]),NA),ifelse(sum(sel)>0,max(transport$transport.nt3[sel]),NA)))
}))
transport$tranport.ngf.longer=tempax[,1]
transport$tranport.nt3.longer=tempax[,2]
transport$diff_transport_longer=tempax[,1]-tempax[,2]
#write.csv(transport,""./zenodo/transport/diff_localisationscore/transport.csv")Let's now check that the axonal localisation score does not depend on the read count in cell body nor on the transcript length.
transport <- read.csv("./zenodo/transport/diff_localisationscore/transport.csv")
transport <- read.csv("./data/transport.csv")
PlotWithAVG <- function(myxval=transport$cb.ngf,myyval=transport$transport.ngf,myxlab="avg.ngf,cb [log2]",myylab="transport efficiency",myylim=c(-4,4),myxlim=c(0,20)){
BY=0.5
plot(myxval,myyval,pch=19,col=rgb(0,0,0,0.1),cex=0.2,xlab="",ylab="",las=1,main="",frame=F,xlim=myxlim,ylim=myylim,cex.lab=0.7,cex.axis=0.7)
grid()
myY<-tapply(myyval,cut(myxval,breaks=seq(from=min(myxval),to=max(myxval),by=BY),include.lowest=T),FUN=function(W)return(mean(W,na.rm=TRUE)))
myX<- seq(from=min(myxval),to=max(myxval),by=BY)[-1]
lines(myX,myY,col="red")
mtext(side=1,line=2,text=myxlab,cex=0.7)
mtext(side=2,line=3,text=myylab,cex=0.7)
}
layout(matrix(c(1,1,2,2,3,4,5,6,7,8,9,10),ncol=4,nrow=3,byrow = TRUE))
hist(transport$transport.ngf,col=rgb(129/255,164/255,214/255,0.5),breaks=20,las=1,cex.lab=0.7,cex.axis=0.7,xlim=c(-4,4),main="",xlab="localisation score [NGF]")
hist(transport$transport.nt3,col=rgb(175/255,113/255,175/255,1.0),breaks=20,las=1,cex.lab=0.7,cex.axis=0.7,xlim=c(-4,4),main="",xlab="localisation score [NT3]")
PlotWithAVG(myxval=transport$cb.ngf,myyval=transport$RUD.ngf,myxlab="avg.ngf.cb [log2]",myylab="log2(AX/CB)",myylim=c(-10,4),myxlim=c(0,20))
PlotWithAVG(myxval=transport$cb.ngf,myyval=transport$transport.ngf,myxlab="avg.ngf,cb [log2]",myylab="transport efficiency",myylim=c(-4,4),myxlim=c(0,20))
PlotWithAVG(myxval=transport$cb.nt3,myyval=transport$RUD.nt3,myxlab="avg.cb [log2] in NT3",myylab="log2(AX/CB)",myylim=c(-10,4),myxlim=c(0,20))
PlotWithAVG(myxval=transport$cb.nt3,myyval=transport$transport.nt3,myxlab="avg.cb [log2] in NT3",myylab="transport efficiency",myylim=c(-4,4),myxlim=c(0,20))
PlotWithAVG(myxval=log10(transport$txLength),myyval=transport$RUD.ngf,myxlab="tx length [log10]",myylab="log2(AX/CB)",myylim=c(-10,4),myxlim=c(2.5,4))
PlotWithAVG(myxval=log10(transport$txLength),myyval=transport$transport.ngf,myxlab="tx length [log10]",myylab="transport efficiency",myylim=c(-4,4),myxlim=c(2.5,4))
PlotWithAVG(myxval=log10(transport$txLength),myyval=transport$RUD.nt3,myxlab="tx length [log10]",myylab="log2(AX/CB)",myylim=c(-10,4),myxlim=c(2.5,4))
PlotWithAVG(myxval=log10(transport$txLength),myyval=transport$transport.nt3,myxlab="tx length [log10]",myylab="transport efficiency",myylim=c(-4,4),myxlim=c(2.5,4))
4 Differential localisation analysis
We can now identify those 3' UTR isoforms which are over-transported in NGF and NT3.
#A. More transported in NGF axons (avoid selecting candidate of axonal remodeling i.e. transport longer should as well be in the same direction)
#Should be more transported in NGF samples
sel_diff <- transport$diff.transport>1.0
#Should be transport or over-transported in NGF
sel_NGF <- transport$transport.ngf>=(1.0)
#Should be under-transported in NT3
sel_NT3 <- transport$transport.nt3<=(0.0)
#If longer, then should be in the same range of differential transport otherwise risk to be axonal remodelling candidate
sel_rem <- is.na(transport$tranport.ngf.longer)#|transport$diff_transport_longer>1.0
#selection for coverage
sel_expression <- transport$anno_tot.NGF.axon.1.raw.corrected.match.rownames.myOut.ngf...>=40&transport$anno_tot.NGF.axon.2.raw.corrected.match.rownames.myOut.ngf...>=40
more_transported_NGF <- data.frame(coord=transport$coord,transport)[sel_diff&sel_NGF&sel_NT3&sel_rem,]#482
#B. More transported in NT3 axons
#Should be more transported in NT3 samples
sel_diff <- transport$diff.transport<=(-1.0)
#Should be transport or over-transported in NT3
sel_NGF <- transport$transport.nt3>=(1.0)
#Should be under-transported in NT3
sel_NT3 <- transport$transport.ngf<=(0.0)
#If longer, then should be in the same range of differential transport otherwise risk to be axonal remodelling candidate
sel_rem <- is.na(transport$tranport.nt3.longer)#|transport$diff_transport_longer>1.0
#selection for coverage
sel_expression <- transport$anno_tot.NT3.axon.1.raw.corrected.match.rownames.myOut.ngf...>=40&transport$anno_tot.NT3.axon.2.raw.corrected.match.rownames.myOut.ngf...>=40
more_transported_NT3 <- data.frame(coord=transport$coord,transport)[sel_diff&sel_NGF&sel_NT3&sel_rem,]
oi<- setdiff(more_transported_NGF$uniqueIDs,more_transp_NGF$uniqueIDs)
View(transport[match(oi,transport$uniqueIDs),])#LS
more_transported_NGF <- read.csv("./zenodo/transport/diff_localisationscore/more_transport_NGF.csv")
more_transported_NT3 <- read.csv("./zenodo/transport/diff_localisationscore/more_transport_NT3.csv")
#Relative Abundance Ratios
more_NGF <- read.csv("./zenodo/transport/diff_transport_ratios/more_in_NGF.csv")
more_NT3 <- read.csv("./zenodo/transport/diff_transport_ratios/more_in_NT3.csv")
par(mfrow=c(2,2),mar=c(4,4,4,4))
colsNGF <- c("#81A4D6","#2D598E","#083872")
colsNT3 <- c("#AE73B1","#79387C","#57055B")
mycols <- rep(rgb(0,0,0,0.1),nrow(transport))
mycols[transport$uniqueIDs%in%more_transported_NT3$uniqueIDs]<- colsNT3[1]
mycols[transport$uniqueIDs%in%more_transported_NGF$uniqueIDs]<- colsNGF[1]
plot(avg.cb[match(transport$uniqueIDs,names(avg.cb))],transport$diff.transport,cex=0.3,pch=19,col=mycols,las=1,xlab="",ylab="",frame=FALSE,xlim=c(3,16),ylim=c(-6,6),main="Analysis based on Localisation Score")
mtext(side=1,line=2,text="average expression CB [log2]",cex=0.7)
mtext(side=2,line=2,text="diff transport efficiency",cex=0.7)
grid()
text(x=10,y=5,col=colsNGF[1],lab=paste("#over in NGF=",nrow(more_transported_NGF)),cex=0.7)
text(x=10,y=-5,col=colsNT3[1],lab=paste("#over in NT3=",nrow(more_transported_NT3)),cex=0.7)
plot(avg.cb[match(transport$uniqueIDs,names(avg.cb))],myOut$dRUD[match(transport$uniqueIDs,myOut$uniqueID)],cex=0.3,pch=19,col=mycols,las=1,xlab="",ylab="",frame=FALSE,xlim=c(3,16),ylim=c(-6,6),main="Analysis based on Localisation Score")
mtext(side=1,line=2,text="average expression CB [log2]",cex=0.7)
mtext(side=2,line=2,text="diff abundance ratios",cex=0.7)
grid()
mycols <- rep(rgb(0,0,0,0.1),nrow(transport))
mycols[transport$uniqueIDs%in%more_NT3$uniqueID]<- colsNT3[1]
mycols[transport$uniqueIDs%in%more_NGF$uniqueID]<- colsNGF[1]
plot(avg.cb[match(transport$uniqueIDs,names(avg.cb))],transport$diff.transport,cex=0.3,pch=19,col=mycols,las=1,xlab="",ylab="",frame=FALSE,xlim=c(3,16),ylim=c(-6,6),main="Analysis based on abundance ratios")
mtext(side=1,line=2,text="average expression CB [log2]",cex=0.7)
mtext(side=2,line=2,text="diff transport efficiency",cex=0.7)
grid()
text(x=10,y=5,col=colsNGF[1],lab=paste("#over in NGF=",nrow(more_NGF)),cex=0.7)
text(x=10,y=-5,col=colsNT3[1],lab=paste("#over in NT3=",nrow(more_NT3)),cex=0.7)
plot(avg.cb[match(transport$uniqueIDs,names(avg.cb))],myOut$dRUD[match(transport$uniqueIDs,myOut$uniqueID)],cex=0.3,pch=19,col=mycols,las=1,xlab="",ylab="",frame=FALSE,xlim=c(3,16),ylim=c(-6,6),main="Analysis based on abundance ratiose")
mtext(side=1,line=2,text="average expression CB [log2]",cex=0.7)
mtext(side=2,line=2,text="diff abundance ratios",cex=0.7)
grid()